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Abstract 

The lattice Bhatnagar-Gross-Krook (LBGK) model has become the most popular one in the lattice Boltz¬ 
mann method for simulating the convection heat transfer in porous media. However, the LBGK model 
generally suffers from numerical instability at low fluid viscosities and effective thermal diffusivities. In 
this paper, a modified LBGK model is developed for incompressible thermal flows in porous media at the 
representative elementary volume scale, in which the shear rate and temperature gradient are incorporated 
into the equilibrium distribution functions. With two additional parameters, the relaxation times in the 
collision process can be fixed at a proper value invariable to the viscosity and the effective thermal diffu- 
sivity. In addition, by constructing a modified equilibrium distribution function and a source term in the 
evolution equation of temperature field, the present model can recover the macroscopic equations correctly 
through the Chapman-Enskog analysis, which is another key point different from previous LBGK models. 
Several benchmark problems are simulated to validate the present model with the proposed local computing 
scheme for the shear rate and temperature gradient, and the numerical results agree well with analytical 
solutions and/or those well-documented data in previous studies. It is also shown that the present model 
and the computational schemes for the gradient operators have a second-order accuracy in space, and better 
numerical stability of the present modified LBGK model than previous LBGK models is demonstrated. 
Keywords: 

Lattice Bhatnagar-Gross-Krook, numerical stability, Generalized model, Convection heat transfer, Porous 
media 


1. Introduction 

Fluid flow and convection heat transfer in porous media have long been a subject of research due to its 
comprehensive relevance in engineering and scientific applications, such as geothermal energy extraction, 
heat exchangers and electronic cooling instruments, chemical catalytic reactors, and pollution transport 
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and mineral processing. Over the past several decades, considerable investigations and applications have 
been devoted to the convection heat transfer in porous media by many researchers using various traditional 
numerical methods (eg. finite difference, element and volume methods). A comprehensive review on this 


subject has been made by Cheng [l|, Nield and Bejan 

As a powerful computational tool based on the kinetic theory, the lattice Boltzmann method (LBM) 
has been successfully applied to simulating complex fluid flows and modeling the physics in fluids 

Some attractive advantages of LBM over the traditional numerical methods have been explained in 


Q| and Vafai Q . 


references [t{| and [§|. Thanks to the mesoscopic and kinetic nature, the LBM has been widely applied to 


. 


the fluid flow and thermal problems in porous media la tip, ii] after its emergence |l2(. Generally speaking, 
applications of the LBM to porous flows in the literature center around two scales: the pore scale and the 
representative elementary volume (REV) scale. At the pore scale [l2, Hi], 14, lIITgI. [17I] . the standard lattice 
Boltzmann equation (LBE) is used to simulate fluid flows in pores, and the local information of flow can 
be directly obtained. Therefore, the LBM at the pore scale can be severed as the most straightforward way 
to investigate the macroscopic relations and reveal the microscopic mechanism of porous flows. However, 
the detailed geometric information of the pores is needed for this approach, and thus the computational 
domain size cannot be too large in view of the limited computer resources. An alternative approach is to 
investigate the averaged quantities at the REV scale. In the LBM at the REV scale, an additional term 
based on some semi-empirical models is incorporated into the standard LBE to account for the presence of a 
porous medium .9, [10|, [llj. Evidence from the literature has demonstrated the REV approach to be simple 
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2lJ. 


and computationally efficient for modeling flows and transport problems in porous media [18 
Guo and Zhao Q] proposed a lattice Boltzmann (LB) model for simulating fluid flow in porous media, 
in which the porosity is introduced in the equilibrium distribution function (EDF), and a forcing term is 
included in the LBE to account for the linear and nonlinear drag forces of the porous media. Further, Guo 


and Zhao flfl ] extended the isothermal LBM to thermal flows in porous media by adding a temperature 
distribution function for the evolution of temperature field. Subsequently, Seta [11] confirmed the reliability 
and computational efficiency of the LBM in simulating natural convection in porous media. Shokouhmand 


et al. 


22| conducted simulations on laminar flow and convective heat transfer in conduits partially and fully 


filled with porous media. Rong et al. H proposed a LB model particularly for axisymmetric thermal flows 


in porous media. Abrach et al. 


24] employed the LBM for the heat and mass transfer during drying of 


deformable saturated porous media. Recently, Gao et al. [25| developed a thermal LB model for simulating 
the non-equilibrium natural convection problems in porous media. 

All the mentioned-above LB models for porous flow and heat transfer problems at the REV scale are 
based on the Bhatnagar-Gross-Krook (BGK) collision model [3- Among these models, both the viscosity 


of fluid and the effective thermal diffusivity of solid matrix are directly determined by the relaxation times. 
The most well-known criticism on these BGK-based models is the numerical instability for moderately low 
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viscosity and/or effective thermal diffusivity. One alternative option to resolve the shortcomings of the BGK 
model is to employ the multiple-relaxation-time (MRT) collision model 27, 28, 3|. Very recently, such effort 
has been made by Liu et al. H who constructed a MRT-LB model within the double-distribution-function 


(DDF) framework. They showed that numerical stability of the BGK model is well improved by the MRT 
model when simulating convection heat transfer in porous media. On the other hand, compared with the 
MRT model, the BGK counterpart has become the most used form of the LBM owing to its extreme 
simplicity and computational efficiency. Therefore, it is of much significance for us to develop a new lattice 
BGK (LBGK) model which can enhance numerical stability at low viscosities and thermal diffusivities. 


Within the framework of LBGK model, Inamuro m proposed a lattice kinetic scheme (LKS) for simu¬ 
lating fluid flows with heat transfer in the absence of porous medium. In this method, an additional term 
relating with the shear rate (and temperature) is incorporated into the EDF to represent the fluid viscos¬ 
ity (and the thermal diffusivity). The relaxation time is fixed at unity while the fluid viscosity (and the 
thermal diffusivity) is determined by another parameter independent of the relaxation time. Thus, better 
numerical stability than the standard LBGK model can be achieved by the LKS at relatively small viscosity 
(and thermal diffusivity). Later, this LKS was extended to develop LBGK models for two-phase fluid flows 




MIM and non-Newtonian fluid flows [jyj]. However, in these LKS schemes, the mass conservation is not 


satisfied unless the fluid density is constant [ 35 I > and the gradient operators of velocity and temperature oc¬ 
curred in the EDF are calculated by a finite-difference scheme. This non-local treatment not only spoils the 
computational efficiency, but also brings about some difficulties to the implementation of complex boundary 
conditions with a local lattice scheme. Although Peng et al. [36] computed the stress tensor locally in vis¬ 
cous thermal flows, the problem of mass nonconservation still remains. Recently, Wang et al. H improved 
the original LKS to develop a mass-conserving and localized LBGK model for non-Newtonian fluid flows. 
However, as far as we know, no works have been reported on extending the LKS to improve the numerical 
stability of LBGK models for incompressible thermal flows in fluid-saturated porous media. 

In this work, we propose a modified LBGK model for simulating convection heat transfer in porous media 
at the REV scale. By introducing the shear rate and temperature gradient into the EDFs as the original 
LKS, the dimensionless relaxation times can be fixed at proper values invariant to the fluid viscosity and 
thermal diffusivity which are determined by two additional parameters. In addition, with a modified EDF 
and a source term in the evolution equation for the temperature field, the macroscopic equations can be 
correctly recovered from the present LBGK model through the Chapman-Enskog analysis. In the following, 
the modified LBGK model is first presented, and the Chapman-Enskog analysis is then given to recover the 
macroscopic equations. Subsequently, a local scheme, instead of the non-local finite-difference schemes, is 
presented for computing the shear rate and temperature gradient. Finally, some benchmark numerical tests 
are carried out to validate the present model. The numerical results show that the present modified LBGK 
model and the local scheme for the gradient operators are both second-order accurate in space. It is also 
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confirmed that the present LBGK model are much more stable compared with the standard LBGK model. 


2. Macroscopic equations 


For a homogeneous, isotropic and fluid-saturated porous medium, the local thermal equilibrium assump¬ 
tion between the fluid and solid is invoked Q. The viscous fluid flow in porous media is assumed to be 
incompressible and the Boussinesq approximation is valid. Neglecting viscous heat dissipation and com¬ 
pression work done by the pressure, the macroscopic gov erning equations for the convection heat transfer in 
porous media at the REV scale can be written as 
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30 


37| 


V • u = 0, 

757 + (u • V) = —-V(ep) + i/ e V 2 u + . 

ut \ £ / On 


dT 


a— + u ■ VT = V • (a e VT) + Q, 


(la) 

(lb) 

(lc) 


where po is the mean fluid density, w, p and T are the volume-averaged fluid velocity, pressure and tem¬ 
perature, respectively; £ is the porosity of the porous medium, v e is the effective kinematic viscosity, and 
Q is the internal heat source term; a = [ epfC p f + (1 — £)p s c ps \/(pfC p f ) is the thermal capacity ratio be¬ 
tween the solid and fluid phases, with pf(p s ) and c p f(c ps ) being the density and specific heat of fluid (solid) 
phase, respectively; a e is the effective thermal diffusivity, and relates with the effective conductivity (fc e ) as 
a e = k e /(pfC p f). F represents the total body force stemming from the presence of a porous medium and 
other external force, and are expressed as 


F =-^u-^j=\u\ u + £ G , ( 2 ) 

where v is the kinematic viscosity (not necessarily identical to v e ), K is the permeability, F e is the geometric 
function, and G denotes the external body force. With the Boussinesq approximation, the body force G 
encompassing the buoyancy force and other external forces is described by 


G = g[3(T - T 0 )j + a , 


(3) 


where g is the gravity acceleration, /3 is the coefficient of thermal expansion, Tq is the reference temperature, 
j is the unit vector in a direction opposite to gravity, and a is the acceleration induced by other external 
force. Based on Ergun’s empirical relation, F e and K can be written as 

F=-^~ K - ^ (4) 

- v'T 50?’ 150(1-e) 2 ’ l ' 

where d p is the diameter of filling solid particle. 

The convection heat transfer problems governed by Eq. © can be characterized by several dimensionless 
parameters: the Darcy number Da , the viscosity ratio J e , the Reynolds number Re (for mixed convection 
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flow), the Rayleigh number Ra (for natural convection), the internal Rayleigh number Raj (for thermal 
convection flows with internal heat source), the Prandtl number Pr, which are defined as follows: 


^ K 
Da =L ?’ 


J — — 

V 


Re = 


LU 


Ra = 


g/3ATL 3 


Rai = 


gPQL 5 


Pr = —, 

O/.Q 


( 5 ) 


v va e va* 

where L is the characteristic length, U is the characteristic velocity and AT is the characteristic temperature 
difference. 


3. Modified LBGK model for thermal flows in porous media 

In this section, we follow the DDF method to develop a thermal LB model for fluid flow and convection 
heat transfer in porous media at the REV scale. The model is based on the BGK collision operator, and 
is constructed with the idea of the LKS. In this model, the flow field is modeled by a LBE of the density 
distribution function, and the temperature field is modeled by another evolution equation of the temperature 
distribution function. For the shear rate and temperature gradient, a local computational scheme is presented 
especially. 


3.1. Lattice Boltzmann equation for the flow field 

For the pressure and velocity fields governed by Eqs. (flail and (00), the evolution equation of the modified 
LBGK model is given by 


fflx + dSt,t + S t ) - fflx,t) = —— \fi{x,t) - f^ eq \x,t)] + 5 t Fflx,t), 

Tf L J 


( 6 ) 


where fflxfl) is the density distribution function for the particle with velocity d at time t and position x, 
St is the time increment, Tf is the dimensionless relaxation time, f^ eq \x,t) is the equilibrium distribution 
function, and F t is the external force term. 

In this work, the original LKS for fluid flows in plain media is extended to incompressible fluid flows in 
porous media, in which the above EDF is defined as [fioj, 31. 35. 


/!”>(*,«) = 


Po - (1 - w 0 )§£ + PoSo(u) + por 0 (u), i = 0 
+ Posflu) +p 0 n(u), ifl 0 

where io % is the weight coefficient, c s is the speed of sound, and 


i(u) = i 


d ■ u uu : ( CjCj — c 2 s I ) 


2 eci 


flu) = 


AS t S : (dd - c 2 s I) 

1 2c? 


( 7 ) 


( 8 ) 


Here, I denotes the identity tensor, and the shear rate S (S = Vtt + (Vu) ) together with an additional 
parameter A are included in the EDF. For the forcing term, the appropriate form of T) to achieve correct 


hydrodynamic equations is taken as 


lOUUl: 


Ft = iOipa ( 1 - — 


Ci-F (uF + Fu) : (cjC* - c 2 I) 


2 eci 


( 9 ) 
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The pressure and fluid velocity are defined as 


P = 


e(l - w 0 ) 


^2 fi + Tf5 t F 0 + p 0 s 0 {u) , 

iji o 


u =—y Cifi + — f 
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( 10 ) 


Due to the nonlinear relation of F with u in Eq. ©, the fluid velocity u is calculated by a temporal velocity 
v and is given by 


v 1 

u = - =. v = — 

Co + V c o + c iM Po 

where the two parameters cq and ci are respectively given by 


y c ifi + -^g, 


( 11 ) 



1 + e 


St”_ 
2 K 


Cl 


St F e 
£ 2 \[K 


( 12 ) 


Note that implementation of the above LBGK model depends on the underlying lattice. For conve¬ 
nience of study, the flow problems considered in this paper are limited to be two dimensional. But it is 
straightforward for us to extend the present model to three dimensional case. For the D2Q9 model, the 
discrete nine velocities are given by c 0 = 0, Cj = c(cos((*— 1)7t/2), sin((i — l)7r/2)) for i = 1 — 4, and 
a = V2c (cos ((* — 5)7t/2 + 7t/4), sin((* — 5)n/2 + 7t/4)) for i = 5 — 8. Here, c = 5 x /5 t is the lattice speed 
with S x denoting the lattice spacing. Accordingly, the sound speed c s = c/y/ 3, and the weight coefficients 
are given by wo = 4/9, Wj = 1/9 for i = 1 — 4, and Wj = 1/36 for % = 5 — 8. 

Through the Taylor expansion and the Chapman-Enskog analysis, the macroscopic equations (Hal) and 
CB can be derived from the LBE (0 (The details are presented in | Appendix A| ), and the effective kinematic 
viscosity is determined by v e = c 2 s (rf — A— 1/2 )6t- We would like to point out that when A in the EDF (j8j) 
is set to be zero, the present LBE © is simplified to the standard LBGK model in which the viscosity v e 
is directly determined by Tf, i.e., u e = c 2 (rf — 1/2 )5t- 


3.2. Lattice Boltzmann equation for the temperature field 

The governing equation C3 is a convection-diffusion equation (CDE) with a heat source term, in which 
the velocity u obeys the incompressible generalized Navier-Stokes equations Gil) and m- For most of 
the existing LBGK models for CDE Gib several unwanted deviation terms are ignored with the additional 
condition, to L/c s (to and L are the characteristic time and length, respectively) in the Chapman-Enskog 
analysis to derive Eq. (flcl) (cf. Ref. [k]] and references therein). In order to eliminate this imperfectness, as 


well as inspired by the idea of Chai and Zhao 


, a modified EDF and a source term should be constructed 


in the evolution equation of the temperature field. Additionally, it is noted that that the temperature 
gradient is included in the LKS. Thus, the evolution equation of the present modified LBGK model for CDE 
GB is written as follows 

gi(x + Cidt,t + S t ) - gi(x,t) = gi(x,t) - g\ eq) (x,t) + 6 t Pi(x,t) + S t Qi(x,t), (13) 

tt L J 
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where tt is the dimensionless relaxation time, gi is the temperature distribution function, and g( q) is the 
temperature EDF and defined by 

Ci ■ u uu : (cjCj — c 2 s I) 


(eg) rr, 

9i = WiT 


2e ci 


+ WiT 


ep 

c 2 s Po 


iBStCi ■ V'/'. 


(14) 


where the coefficient Wi is given by Wi = u>i [i ^ 0), tuq = — The two source terms Pi and Qi are 

taken as 


Pi = OJi |^1 - 

Qi — i ( 1 


1 

2tt 

1 

2tt 


Ci - (TF + epVT/po) 


1 1 ^ uu : (c, VT) 

e er 


1 + 


Q- 


(! 5 ) 

(16) 


It is noted that the discrete source term Qi corresponds to the source term Q resided in Eq. m- 

The LBE (fl3l) is implemented separately with two substeps, i.e., the collision step and streaming step: 

Collision: g*(x,t) = gi(x,t) - \gi(x,t) - g[ eq \x,t)\ + 5 t Pi(x, t) + 6 t Qi(x, t), (17) 

t t L J 


Streaming : g t (x + CiS t ,t + S t ) = g* ( x , t), 


(18) 


where g* (x, t) denotes the postcollision distribution function. At each time step when the two substeps are 
completed, the macroscopic temperature T are computed by 

St. 


°T = '^g l 


-Q. 


(19) 


In the LBE (TTdl) for the temperature filed, both the EDF gf q and the source term Pi depend on e and a with 
which the influence of porous media is introduced. It should be noted that although a linear EDF is widely 
employed in previous LBGK models for Eq. m 0 Sj Si ) large errors will be encountered in solving 
the CDE (fTcl) by these models when the resulted numerical diffusion coefficient from the Chapman-Enskog 


analysis is proportional to the velocity square 


39 


41]. Additionally, when cr = e= l,B = 0 and without 


the forcing term Qi, the present LBGK model will reduce to that proposed by Chai and Zhao 


39] for the 


CDE without porous media, which can thus be considered as a special case of the present model. 

By the Chapman-Enskog analysis on the LBE (1131) . the temperature equation © can be recovered 
exactly. In what follows, the detailed mathematic derivations will be shown. To this end, the Chapman- 
Enskog expansions are first applied to the distribution function gi , the derivatives of time and space, and 
the internal heat source term Q: 

n - „ (0) 4- \n W + \ 2 „ (2) J- 

9i — 9i + A9i + A 9i H-> 


dt — A dt 1 + A 2 <9 t 2 i 
Q = AQ (1) , 


V = AVi 
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(20a) 

(20b) 

(20c) 






















where A is a small expansion parameter having the magnitude of the Knudsen number. Since the spatial 
derivative of temperature is contained in the EDF g^ eq ' > and the force term Pj, the following multiscaling 


expansions are also introduced 


Q: 


g^=gt W + Xgt W , P, = 


where g^°\ 9i^\ and P can be explicitly written as 

Cj • u uu : (ciCi — c 2 s I) 


e(0) rri 

9i = UiT 


2 ecj 


+ w z T 


£ P 

~~0 5 

cjpo 


gf 1] = aUiBStd ■ Vi T, 


P (1) = Ui I 1 - 


2tt ) 


Ci ■ ( TF + £pViT/po) , (1 1\ uu : (c,ViT) 

e a 


( 21 ) 

(22a) 

(22b) 

(22c) 


After some computations, we can get the following velocity moments 


E e(0) rp e (°) rp e (°) Tuu 2 rp T . £ P T T e(0) 2rp A 

9i K, =aT, l^dgC ’ =Tu, 2^ c i c ^9i = -h c s aTI-\ - 1, 2^dc i dg l = c s TA ■ u 


P o 


5Z 9i =0, Ci 9i = (rc 2 s B6 t S7 iT, ^ dd9^’ = 0, 


e(l) 


Ed‘>=o. A>rf> = (i-J- 


TF 


£pViT /1 1 


Po 


- uu ■Vi T 

£ (7 




2rr / a 


(23a) 

(23b) 

, ^ c,c^ (1) = 0, (23c) 

i 

(23d) 


where A • u = u a 5p-y + Up5 ai + 3 , in which 5 a p denotes the Kronecker delta with two indices, and 

Q[ 1] = A Qi. 

The Taylor series expansion applied to Eq. m up to second order in S t leads to 


S t 2 1 

Di9i + —D.gt = [9i - 9i 


(eq) 


+ $tPi + S t Qi, 


(24) 


where Di — d t + d • V. Substituting Eqs. and (HHll into Eq. (Ell . one can obtain the consecutive orders 
of Eq. (l24l) in terms of A: 


\0 ■ „ (0) - „ e(0) 
A ■ 9i ~ 9i > 

\1 . n JO) _ 

A ■ u li9i — — 


TrS t - 


(1) e(l) 

9i ~9i 


+ PV+QW, 


\2 . f) JO) , n J 1 ) , n 2 « (0) - 1 « (2) 

A ■ OtaPi +^ L> li9i — ) 

where Ph = 9 tl + c,; • Vj. With Eq. (I251)f) . Eq. (I25cl) can be rewritten as 

fl JO) , A 1 \ n Jl) , 1 n e(l) <5t n /p(l) , n(lA_ 1 J2) 

^ + ^ D u9i +2 Du \ Pi +Q * )-~^5 t 9i ‘ 


(25a) 

(25b) 

(25c) 


( 26 ) 
























Note that aT = ]T\ gf q) = JA ^ + 4fQ. Based on Eqs. (l2l1) - (E3ll . it is easy from Eq. (I25al) to obtain that 

£ a (,> = -£ei 1> , EA-o^n. ( 27 ) 

i i i 

Along with these derived equations, by taking summations of Eqs. (I25bl) and (l2fil) over i, the macroscopic 
equations at the fi = At and f 2 = A 2 t time scales can be obtained 


d tl {oT) + V 1 -(Tu) = QW, 


(28a) 


dtJvT) + (l - J-'j Vl • ^c^'j + • (ac 2 s B6 t V 1 T) 


+ 1 - 


2t t / 2 


-Vi • TF (1 > + 


epWiT (1 1 


P o 


e a 


uu • ViT H- Q 

<7 


( 1 ) 


= 0. 
(28b) 


From Eq. (I25bl) . g can be draw out to further compute )TA cig^\ After some tedious but standard 
algebraic manipulations, we can get 

iT N 


- -*[^) + V,. (^) + V, (f) - (1 - £) (-<« + ^ 


* 1 2t t ) (e cr * UU ' VlT * 1 




Q 


(i) 


(29) 


+ c 2 s a (B - t t ) 6tVxT. 


2 tt ) <J 

With the aid of the incompressible Navier-Stokes equations at t\ time scale (see|Appcndix A|) 


Vi • u = 0, 


d tl u + Vi ■ ( —) = —-Vr(ep) + F (1) , 


V e ) 


P o 


(30a) 

(30b) 


and Eq. II 28 all , the first two terms in the square bracket on the right hand side of Eq. (E 2 D can be evaluated. 
After some standard algebra, it can give us 


a tl (T«) + Vi 


<cr) 


d u 

(Tu) + Vi 

/ Tuu 

) + 


V 

J 



/uu\ 

u 

T 

d tl u + Vi 

■ — 

+ - 



V £ J\ 

a 

tfM TVl 

(cp) + u Q m 


1 1 

e a 


Vi • ( Tuu ) 


uu ViT. 


1 1 

£ a 


ViT (31) 


Further, application of Eq. m to Eq. (l29l) yields 


E c ^ (1) = - 


TF [ 


(1) , £ P VjT \ 5 t (1 1 


P o 


St u . 


uu • VlT - + cia (B - t t ) S^iT. (32) 

2 <7 


As a consequence, Eq. (I28bl) can be rewritten as 


dt 2 (crT) + Vi • 


aci[B-T T + -} StViT 


= 0. 


(33) 


9 



























Combining the derived equations at t\ and £2 time scales (Eqs. (1 28 al) and ( 1551 ) 1 . the final CDE can be 
recovered 

d t (aT) + V • (Tu) = V • (a e VT) + Q, (34) 


where a e is the effective thermal diffusivity and determined as 





(35) 


Under incompressible flows condition, i.e., V • u = 0, it should be seen that the derived CDE ( 1551 ) is just 
Eq. CEB. 

From the presented Chapman-Enskog analysis, it is clearly seen that the macroscopic equations Eq. © 
for convection heat transfer in porous media can be exactly recovered from the present modified LBGK 
model. Now, some comments on the present modified LBGK model are in order. First, the fluid viscosity 
and the thermal diffusivity are not determined only by the relaxation times. Thus, as compared with the non- 
modified LBGK model, the present LBGK model can achieve better numerical stability and boarder range 
of Prandtl number values in the simulations. Second, the EDF g^ q (Eq. ( 1151) 1 in the present LBGK model is 
nonlinear in terms of u , and thus avoids the numerical diffusion coefficient resulted in the recovered CDE. 
Third, the temperature gradient VT appearing in g[ eq ^ and Pi can be computed by a local scheme instead 
of a finite difference scheme (as will be presented in the forthcoming section). Thus, the locality of collision 
process (Eq. ( 1571) 1 of the present model is preserved. Finally, as e — > 1 , one can see that the macroscopic 
equation Eq. m will be that for free convection heat transfer without porous media. Accordingly, the 
present LBGK model is actually reduced to that for thermal convective flows in the absence of porous 
media. 


3.3. Local scheme for the shear rate and the temperature gradient 

In the framework of LBM, it is direct to use the traditionally nonlocal finite-difference schemes to compute 
spatial gradients. However, the shear rate and the tem per ature gradient can be locally calculated by the 


nonequilibrium part of the distribution function 


39 


43 


44j without the influence of porous media. In this 


work, we will provide a local computational scheme for these two gradient terms in the presence of porous 
media, which has not been reported in the literature yet. 

First, the local computing scheme for the shear rate is derived. Reviewing Eq. (IA.13I) in |Appendix A 


Y CidfY = c 2 s po(A - Tf^tSx - ^p 0 


S t (uFA) fA) 


u 




and multiplying £ on its both sides with the relations ffY = /,; — f- 0 ^ + 0(£ 2 ) and f- 0 ^ = f^°\ we can 
obtain that 


S = 


Si C i C i 


fi-ft 


( 0 ) 


I St n l uF 1 Fu \ 

+ T Po {— + —) 


c 2 s p 0 (A - T f ) S t 
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(36) 



















where is given by Eq. 0. and its velocity moments can be found in Eq. (IA.4I) . 

Next, we turn to the computational scheme for the temperature gradient VT. Similar to the shear rate, 
multiplying Eq. ( 1551 ) with A on both sides and recognizing that A g^ = gt — + 0{ A 2 ) and g^ = g^°\ 

we can get 




9i ~ 9, 


e(0) 


f( TF W«) 


ep5 t 

2po 


VT “ 77 ( - “ “ ) uu ' VT + & ( B - t t) StVT. (37) 


Due to the tensor uu contained in the above equation, we cannot directly obtain VT. Fortunately, based 
on the equality VT = I ■ VT, Eq. ( 1571 ) can be put into the following tensor form 


epSt T S t f 1 ll 2 tT3 \ X T 


-VT = ■£' 


9i - 9; 


e(0) 


+ I( TF + ? Q )- (38) 


Now we can obtain VT by solving Eq. ( 1551 ) . Mathematically, this can be done by solving the following 
equation in the matrix form 

(M)(VT) = (N), (39) 


where (M), (VT) and (N) are three matrixes, which are respectively formed by the elements of tensors 


-Tt 1 ~ IQ - t) «« + & (B - t t ) S t I , VT and c, 


/ M\ £pSt A St ( 1 1 


9i ~ 9 


"( 0 ) 


f (TF + hq) 


yU /3 + c 2 s a {B - t t ) S t 5 a p, 


(X7T)p = dpT, 
(N)a = y~]cj, 


9i - 9, 


e(0) 


+ l( TF " + »- 


By some simple computations, the determinant of the matrix (M) can be obtained 

epS t 


det((M)) = 


c;a (B - tt) S t - 


2po 


epSt 

2po 


£ 

2 


(1 

1 \ 

I i2 

— 


1 u 1 


°7 



(40) 

(41) 

(42) 

(43) 


From Eq. ( 1551 ) . one can see that B < tt- Furthermore, in the simulations of this paper, er is set to be 1, 
while £ is varied within 0 < e < 1. Hence, 1/e — 1/er > 0. With these inequalities, we can conclude from Eq. 
( 1551 ) that det((M)) > 0, that is, (M) is an invertible linear matrix. Consequently, (VT) can be determined 
by (VT) = (M)~ 1 (N). 

It is clear that the above formulae to compute the shear rate and the temperature gradient are local 
schemes while no finite difference schemes are employed. This ensures that the collision process of the 
present LBGK model can be locally performed. A few remarks regarding the local schemes for computing 
the shear rate and the temperature gradient are given here. First of all, the local computing schemes are 
still available for the case without porous media only by taking the porosity e = 1. Second, as A and B are 
equal to zero, the shear rate S disappears from the present model, however, the temperature gradient VT 
still needs to be computed since it remains in the source term Pi. Third, when e = er or u • VT = 0 (this 
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corresponds to the pure diffusion system), the second term on the right hand side of Eq. (1371) equals to zero. 
As such, the temperature gradient can be directly obtained, and is given by 


VT = 


Ei c "' 

e(0) 

9i~9i 

+ |(TF+hq) 


( B ~ -t) - ^ 

<5* 


(44) 


3-4- Boundary conditions 

In this work, to specify physical boundary conditions for the distribution functions, we will extend 
the nonequilibrium extrapolation scheme (NEES) [ 45 I to the present LBGK model due to its simplicity, 
second-order accuracy, capability for different boundary conditions, and good robustness. 

The basic idea of NEES is to decompose the distribution function at a boundary node x b into the 
equilibrium part and the nonequilibrium part. Specifically, for the density distribution function /), the 
decomposition at x b is described as 


fi(x bl t) = rf eq \x b ,t) + f[ neq \x b ,t). 


(45) 


The physical boundary conditions are embodied through the equilibrium part, while the nonequilibrium 
part is approximated by certain extrapolation schemes. With the same idea in the present LBGK model, 
the NEES is first generalized for the velocity boundary condition where u(x b ,t) is known but p(x b ,t) is 
unknown. From Eqs. 0 and (|g|). one can see that the pressure p(x b ,t) and the shear rate S(x b ,t) need 
to be evaluated to determine the equilibrium part. For this purpose, we use the pressure p(xf,t) and the 
shear rate S(xf,t) at Xf to respectively replace p(x b , t) and S(x b ,t), and then the EDF f- eq \x b ,t) at x b 
is approximated as 


fi eq \x b ,t) = 


Po - (1 - uo) ep ^i’ t) + PoSo(u(x b: t)) + p 0 r 0 (u(x f ,t)), i = 0 
+ poSi(u(x b ,t)) + pon(u(x f ,t)), i ^ 0 


(46) 


where Xf is the nearest neighbor fluid node away from x b , and 


Si(u(x b ,t )) = Wi 


Ci ■ u(x b ,t) u(x b ,t)u(x b ,t) : ( CiCi - c 2 s I) 


+ 


2 ec\ 


, ri(u(x f ,t)) =Ui 


AS t S(xf,t) : (cjCj - c 2 s I) 
2 cl 


(47) 

Note that S(xf,t) can be calculated directly according to Eq. (13^1) . For the nonequilibrium part at x b , 
f!> neq \x b ,t) is approximated by the nonequilibrium part of fi(xf,t ) 


f( neq) (x b , t) = fi ( x f ,t)~ fl; eq) ( x f ,t ). 

As a whole, the distribution function at the boundary node x b is calculated by 

fi{x b ,t) = fl eq \x b ,t) + \fi(x f ,t) - rf eq] {x f ,t) . 
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(48) 


(49) 
















Thermal boundary conditions with the NEES can be treated in a similar way. Two cases of temperature 
boundary conditions are taken into account. For the case that the temperature T(x b ,t) at x b is only known, 
determination of VT(x b ,t) is needed for g[ eq \x b ,t) (see Eq. (fl4l) b By approximating VT(x b ,t) with 
S7T(xf 1 t) and combining the velocity boundary condition, the equilibrium part g\ eq \x b: t) is specified as 


g[ eq) = UiT(x b ,t) 


Cj ■ u(x b ,t) u(x b ,t)u(x b ,t) : {cjCj - c 2 s I) 
c 2 s 2ercf 

+wiT(x b ,t) £P ^ f,t ^ + aujiB5 t Ci ■ VT(s/,f). 
CtPO 


(50) 


As such, the temperature distribution function gi(x b ,t) is calculated by 


gi{x b ,t) = g[ eq \x b ,t) + 


gi(x f ,t)~ g\ eq) (x f ,t) 


(51) 


For the other case when the temperature gradient VT( x b , t) is known, the local temperature T(x b , t ) at the 
boundary node x b can be taken from the following second-order finite difference approximation 


T(x b t) — ^( x fA) 2A ■ \7T(x bl t) 


(52) 


where Xff is the nearest fluid node of Xf such that A = Xf — x b = Xff — Xf. Hence, the EDF g\ eq \x b ,t) 
at x b can be approximated by 


g[ eq \x b ,t) = uji 


4 T(x f ,t) - T(xff, t) - 2A ■ VT(x b ,t) 
3 


<T + 


Ci ■ u(x b ,t ) 


c 


2 

s 


u(x bl t)u(x b , t) : (CjCj - c 2 s I) 

2 £c s 


+ WiT(x b ,t) £P ^ b,t ^ + atOiBStCi ■ VT(x b ,t), (53) 

c s Po 

and then the distribution function gi(x b ,t ) at x b is calculated according Eq. (l5ll) . Finally, we would like to 
point out that other boundary schemes can also be developed in a similar way for the velocity and thermal 
boundary conditions in the present LBGK model. 


4. Results and discussions 


In this section, we shall validate the proposed LBGK model by some two-dimensional convective heat 
transfer problems in porous media, and compare the simulation results with analytical and previous numer¬ 
ical results. To ensure the numerical results to reach the steady state, the following convergent criterion is 
used in the simulations 


E* II 0(®, t) - ^i. x i * ^ iOOtft) || „ , n ,, , n _ 6 
EJI#M)II <LOxl ° ’ 


(54) 


where (j) is u or T, and its norm is computed respectively with the L 2 and L 1 -norm. Unless otherwise 
specified, values of the following parameters are used: po = 1.0, J e = 1.0, a = 1.0, and a e /af = 1.0 (a/ is 
the thermal diffusivity of the fluid), and the dimensionless relaxation times ry and tt are set to be unity. 
Due to the positivity of the viscosity and effective thermal diffusivity, it should be clear that values of A 
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and B should satisfy the constraint of A < Tf — 1/2 and B < tt — 1/2. As noted for the relaxation times 


elsewhere 


3C 


42|, A and B can also be determined via Pr , Ra and Ma , where Ma is the Mach number 


and defined as Ma = U/c s (U = y/g[3ATL is the characteristic velocity), which are given by 

L l3Pr „ 1 Ma L 


A = T f - Ma 


3 Pr 
~Ra : 


B = tt - 

2 


a 5 X V PrRa 


(55) 


To make the flows fully lie in the incompressible regime, the Mach number should be small and is set to be 
0.1 in this work. 

In addition, the accuracy of the present model and the computing scheme for the shear rate and tem¬ 
perature gradient will also be investigated. Correspondingly, the following relative global error is used 

E* II ®a(x,t) - <$> n {x,t) || 


E{$) = 


(56) 


EJI *»(*,*) || 

where $ represents a scalar or vectorial variable in terms of velocity u and temperature T, and the subscripts 
a and n denote the analytical and numerical results. 


4-1- Mixed convection in porous channel 

The first test problem is the mixed convection flow in a channel filled with a porous medium. The 
schematic of this problem is shown in Fig. [Q The channel of height H and length L is filled with a porous 
media of porosity e. The top wall of the channel is held at a high temperature Th and moves with a uniform 
velocity uq along the x-direction, and the static bottom wall is held at a low temperature T c . A constant 
normal flow of fluid is injected through the bottom wall and is withdrawn at the same rate from the top 
wall. For the injected fluid in the porous channel, it will be sheared by the identical fluid in the normal 
direction of injection. If the nonlinear inertia force (F e = 0) is neglected, the flow at the steady state obeys 


ttttt T = T l 


U 0 


H 


Vo 


Vo T = T r 


ttttt 


L 


Figure 1: Schematic of mixed convective flow in a porous channel. 
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the following governing equations 



u y du x 
e dy 


d 2 


U x £V e 


(57a) 


^% = 9 ^ T ~ T °^ Jc Uy + ay ’ ( 57b ) 

dT 

u y — =V-(a e VT), (57c) 

dy 

where u x and u y are the two components of fluid velocity u , To = (T^ + T c )/2 is the average temperature, 
and a y represents the external force in the y-direction which is given by 


a v R v 0 g(3AT _ 1 


exp (yv 0 /a e ) - 1 


The analytical solutions of u and T in Eq. <53 are identified as 

sinh(C 2 • y/H) 


u x = w 0 exp 


<■(!-) 


T = To + AT 


sinh(C 2 ) 

exp (PrRe ■ y/H) — 1 


Uy - 1 ^ 0 ? 


(58) 


(59a) 


(59b) 


exp (PrRe) — 1 

where Re is the Reynolds number defined by Re = Hvq/v, AT = Th~T c is the temperature difference, and 
the two parameters and £2 are respectively given by 


Re 




2 eJe 


C 2 


1 

2e J e 


Re 2 


4e 3 J e 
Da 


(60) 


From Eqs. II 59 all and (I59bll . the following exact solutions to the gradients of velocity and temperature can 
also be obtained 


du x 

dy 


u 0 


sinh(C 2 ) 


exp 


Mi-0 


^sinh(C 2 • y/H) + ^cosh(C 2 • y/H) 


(61a) 


— = a r PrRe exp {PrRe ■ y/H) 

dy H exp(PrRe) — 1 ' 1 ’ 

In the simulations, the computational domain is chosen to be H x L = 1.0 x 1.0, and some physical 
parameters are listed as follows: e = 0.6, Ra = 100, Pr e = 1.0, Da = 0.01, uq = vq = 0.01. The lattice 
size is N x x N y = 32 x 32, and the values of A and B are determined by A = Tf — 0.5 — Hvo/(c 2 StRe), and 
B = tt — 0.5 — (r/ — 0.5 — A)/(aPrJe) (or B = tt — 0.5 — Hvq /(aRePrc 2 5t)) here. Periodic boundary 
conditions are employed at the entrance and outlet, and the NEES presented in Sec. 13.41 is applied to the 
top and bottom walls for the velocity and temperature boundary conditions. In Fig [2j the comparisons 
of profiles for the velocity and temperature together with their gradients are summarized under different 
Reynolds numbers at Ra = 100, Pr = 1, e = 0.6 and Da = 0.01. As can be observed, the present numerical 
results of velocity and temperature and their gradients agree well with the analytical solutions. 
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Figure 2: Profiles of (a) velocity; (b) temperature; (c) velocity gradient and (d) temperature gradient for different Re at 
e = 0.6, Ra = 100, Pr = 1, Da = 0.01 with a mesh size of 32 x 32. Solid lines: analytical solutions; Symbols: numerical 
results. 


16 
































Figure 3: Relative errors of (a) velocity and (b) temperature fields against lattice spacing 8 X at Re = 5, £ = 0.6, Ra = 
100, Pr = 1, Da = 0.01. The slope of the inserted line is 2.0, which indicates that the present models for the flow and 
temperature fields are of second-order accuracy in space. 


In what follows, the spatial accuracies of the present model for the flow and temperature equations 
are tested by computing the relative errors with different lattice sizes. In Fig. [3l the relative errors in 
velocity and temperature fields are presented at Re = 5, e = 0.6, Ra = 100, Pr = 1, Da = 0.01. The 
results are obtained based on variational lattice spacings from 1/16 to 1/96, and five values of relaxation 
times (t/ = 0.55, 0.8, 1.0, 1 .2, 2.0) are used. As seen from the figure, the present model for the flow and 
temperature fields are second-order accurate in space. In addition, the local computing schemes for the 
shear rate tensor and the temperature gradient are evaluated in terms of spatial accuracy. Note that the 
analytical solutions for the fluid velocity and the temperature have the forms as Eqs. (1 59 all and (I59bl) . Thus, 
the component of (and in S (and VT) is only computed for its relative errors with different lattice 
sizes, which is plotted in Fig. [I] The results shown in Fig. [J] clearly indicate that the present local schemes 
for computing the shear rate tensor and the temperature gradient are of second-order accuracy in space. 

As can be seen from Figs. [3] and 21 the relaxation time has an influence on the computed relative errors 
of velocity and temperature and their spatial gradients, and smaller errors are obtained with the relaxation 
time near unity. This indicates that more accurate results can be obtained when the relaxation time is 
around unity. Next, we turn to investigate the influences of A and B on the spatial accuracy of the present 
model. In the simulations, the relaxation times are fixed at Tf = tt = 1.0, the values of A and B are varied 
in the range of —0.2 < A, B < 0.4, and other parameters are identical to those used before. Fig. [5]shows the 
computed relative errors under different lattice sizes. As indicated in the figure, the variation of A and B 
do not influence the spatial second-order convergence rate of the present model and the computing schemes 
for the shear rate and temperature gradient, but it affects the accuracy of the present model and the local 
scheme. Further, it is found that the influence of A on the accuracy of flow held is different from that of B 
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Figure 4: Relative errors of (a) and (b) against lattice spacing 5 X at Re = 5, e = 0.6, Ra = 100, Pr = 1, Da = 0.01. 
The slope of the inserted line is 2.0, which indicates that the present models for the flow and temperature fields are of second- 
order accuracy in space. 


on the temperature field. For instance, as A = —0.2 and B = —0.2, the relative errors of velocity and shear 
rate are relatively large among the whole computed results (see Figs. [5] (a) and (c)), while the relative errors 
of temperature and its spatial gradient are the smallest (see Figs. [5] (6) and (d)). To reveal the influence 
of these two parameters on the accuracy of computed results more clearly, the four relative errors against 
different values of A and B are plotted in Fig. [HI where S x = 1/32. As can be seen, the relative errors 
of temperature and its gradient increase with increasing B. In contrast, the relative errors of velocity and 
shear rate exhibit two different trends as A increases, that is, they decrease when A is approximately smaller 
than zero, and turn to increase when A is larger than zero. Based on these results, it can be recognized 
that, to derive more accurate results, A and B as well as the relaxation times should be also assigned with 
appropriate values, and especially A should be around zero. 

As noted previously, due to the presence of two additional parameters A and B , the present modified 
LBGK model would be more stable than the standard LBGK model. Now, we will confirm this speculation by 
comparing the results predicted by the present model with those by the LBGK model jn3] at low viscosities. 
In Fig. [7J the computation results with two viscosities v = 9.375 x ICG 6 and v = 9.375 x ICG 7 are shown 


at Re = 5 ,Ra = 100, Pr = l,e = 0.6 ,Da = 0.01, <5^ = 1/32. From Fig. 7(a) one can see that the LBGK 
model y becomes unstable at v = 9.375 x 10 6 (t/ = 0.5009) when t = 46000<5*, while the present LBGK 
model (Tf = 1.0) is stable (see Fig. |7(b)[) and finally yields excellently agreeable results to the analytical 


solution (see Fig. 7(c)). This demonstrates that the present modified LBGK model is more stable than 
the standard LBGK model. To strengthen this statement, we further performed simulations with smaller 




and find that 


values of viscosity, for example v = 9.375 x 10 ' (jf = 0.50009 in the LBGK model 
the present LBGK model can still give excellent agreement results (see Fig. |7(d)[). Additionally, since the 
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Figure 5: Relative errors of (a) and (b) against lattice spacing S x at Re = 5, e = 0.6, Ra =100, Pr = 1, Da = 0.01. 
The slope of the inserted line is 2.0, which indicates that the present models for the flow and temperature fields are of second- 
order accuracy in space. 
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Figure 6: Relative errors of the flow velocity and the shear rate together with the temperature and its gradient at different 
values of A and B. 


Prandtl number is unity in the simulations, we would like to point out that better numerical stability of 
the present LBGK model at low effective thermal diffusivities is actually also verified. These results further 
demonstrate the enhanced numerical stability of the present LBGK model over the standard LBGK model, 
which can be attributed to the independent determination of relaxation times with the viscosity and effective 
thermal diffusivity. 


4-2. Natural convection in a porous cavity 


Natural convection in a square cavity saturated with porous media has been extensively studied as a 
benchmark flow by many researchers [l£J, [ill, |46|, |47], |48|. The schematic of the problem is shown in Fig. [8l 
where the upper and lower walls of the cavity are adiabatic, while the left and right walls are isothermal 
but held at different temperatures Th and T c (7), > T c ), respectively. At all walls of the cavity, the zero 
velocity boundary condition is imposed. The height and width of the cavity are H and L , the temperature 
difference is AT = Th — T c , and the reference temperature is T 0 = (Th + T c )/2. The average Nusselt number 
Nil on the left (or right) vertical wall is defined as 


Nu= J Nu(y)dy , (62) 

where Nu(y) = -~L(dT/dx) wa u/AT is the local Nusselt number. 

We first apply the present modified LBGK model to the case in which e —> 1 and Da tends to infinity 
with 10 3 < Ra < 10 6 . As noted before, the convection heat transfer problem is actually reduced to the 
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u(y,L/2)/u 0 


(a) LBGK [10|, t = 460005*, r f = 0.5009 



(c) Present LBGK, steady-state, ry = 1.0 



(b) Present LBGK, t = 460005*, ry = 1.0 



(d) Present LBGK, steady-state, Ty = 1.0 


Figure 7: Velocity profiles of the mixed convection flow at Re = 5,-Ra = 100, Pr = l,e = 0.6, Da = 0.01 with low viscosities 
v = 9.375 x 10 -6 and v = 9.375 x 10“ 7 . The numerical results are predicted by the present LBGK model and the standard 
LBGK model [lot on a N x x N y = 32 x 32 mesh size. 
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Ux — ti y — 0 



Adiabatic wall 


tt x ti y — 0 


Figure 8: Schematic of natural convection in a fluid-saturated porous cavity. 


case without porous media. In the simulations, we use the grid size of 100 x 100 for Ra = 10 3 , 150 x 150 
for Ra = 10 4 , and 200 x 200 for Ra = 10 5 ,10 6 , respectively. For comparison with previous results, some 
quantities are concerned here: the maximum horizontal velocity component u max at the mid-width (x = Lj 2) 
and its location y max , the maximum vertical velocity component v ma x at the mid-height (y = H/2) and its 
location x max , the maximum Nusselt number Nu max and the corresponding location y^ u , and the average 
Nusselt number Nu along the cold wall. In Table [Q the present results are listed for e = 0.9999, Da = 10 8 , 
Pr = 0.71, and Raj = 0. Here, the computed velocities are normalized by the reference velocity a/L, and 
the x- and y-coordinates are normalized by L and H , respectively. The benchmark data from Refs. |49l |50| 


and the LB results given in Refs. 


30 


are included in the Table for quantitative comparison. It can be 


found that our numerical results are in excellent agreement with those reported data. 

The cases with the influence of porous media are further examined under different Darcy numbers and 
porosities. In Fig. O the streamlines and isotherms predicted by the present model are plotted for the 
Darcy-Rayleigh number Ra* = DaRa = 100 and e = 0.4. In the simulations, the results at Da = 10~ 2 , 
Da = 10 -4 and Da = 10 -6 are obtained with a mesh size of 120 x 120, 200 x 200 and 250 x 250, respectively. 
As can be observed from the figure, for the same Ra* the velocity and thermal boundary layers near the 
hot and cold walls become thinner as Da decreases. At higher values of Da (Da = 10 -2 ), the convective 
mixing is more vigorous inside the cavity, and the isotherms are more sparse near the corners. On the other 
hand, we also find that as Ra increases with a fixed Darcy number, the patters of isotherms change from 
almost vertical to be horizontal in the center of the cavity while being vertical only in the thin boundary 
layers near the hot and cold walls. This indicates that the dominant heat transfer mechanism is varied from 
conduction to convection. All of these observations and findings are in good agreement with those reported 
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(a) Da = 10 —2 , Ra = 10 4 





Figure 9: Streamlines (left) and isotherms (right) of the natural convection in a cavity for e = 0.4, Raj = 0, and Pr = 1.0: 
(a) Da = 10~ 2 , Ra = 10 4 ; (b) Da = 10 ~ 4 ,Ra = 10 6 ; (c) Da = 10 6 , Ra = 10 8 respectively with a mesh size of 120 x 120, 
200 x 200 and 250 x 250. 
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Table 1: Comparisons of the numerical results by the present modified LBGK model with the benchmark solutions I49l. l5(jl and 
the reported LB data Idol , I 51 I (e = 0.9999, Da = 10®, Raj = 0, Pr = 0.71, and grid sizes: 100 X 100 for Ra = 10 3 , 150 X 150 for 
Ra = 10 4 , and 200 X 200 for Ra = 10 s , Ra = 10 6 ). 


Ra(N x x N y ) 


'Umax 

Umax 

Umax 

•Umax 

N Umax 

VNu 

Nu 

10 3 

Ref. [49] 

3.649 

0.813 

3.697 

0.178 

- 

- 

- 

(100 x 100) 

Ref. [51] 

3.6554 

0.8125 

3.6985 

0.1797 

1.5004 

0.90625 

1.1168 


Present 

3.6348 

0.8100 

3.7007 

0.1800 

1.5089 

0.9100 

1.1181 

10 4 

Ref. [50] 

16.1802 

0.8265 

19.6295 

0.1193 

3.5309 

0.8531 

2.2448 

(150 x 150) 

Ref. [30] 

16.1623 

0.8200 

19.6074 

0.1200 

3.5311 

0.8533 

2.2436 


Present 

16.1629 

0.8267 

19.6164 

0.1200 

3.5405 

0.8600 

2.2465 

10 5 

Ref. [50] 

34.7399 

0.8558 

68.6396 

0.0657 

7.7201 

0.9180 

4.5216 

(200 x 200) 

Ref. [30] 

34.7118 

0.8550 

68.4537 

0.0650 

7.7427 

0.9200 

4.5197 


Present 

34.9841 

0.8550 

68.4220 

0.0650 

7.7525 

0.9250 

4.5182 

10 6 

Ref. [50] 

64.8367 

0.8505 

220.461 

0.0390 

17.5360 

0.9608 

8.8251 

(200 x 200) 

Ref. [30] 

65.1052 

0.8500 

219.379 

0.0400 

17.5857 

0.9600 

8.7731 


Present 

65.0494 

0.8500 

218.2964 

0.0400 

17.6373 

0.9650 

8.7792 


Refs 


30 


48 


511 ]. To quantify the results, the average Nusselt numbers of the left vertical wall from the 
present calculations are compared with those previous results, which are listed in Table [2] It is clearly seen 
that the present results agree well with the well-documented numerical results in previous studies. 


4-3. Thermal convection of a heat-generating fluid in a porous cavity with isothermally cooled walls 

In the above tests, the thermal convection problems in a fluid-saturated porous medium is studied but 
without internal heat generation. To further show the potential of the present LBGK model, the natural 
convection flows in a porous cavity with internal heat generation are tested. This kind of thermal convection 
flows have been studied numerically by many researchers 


3(1 . 5*1 , and the convective flows induced by internal 


heat generation in a porous cavity with isothermally cooled walls are considered in this subsection. 

The schematic illustration of simulated problems are shown in Fig. [TO] The walls of the porous cavity 
are maintained at a constant temperature T c and subjected to no-slip boundary conditions. The internal 
heat source term is Q 1 and the temperature difference AT is defined as AT = QL 2 /a e , by which g/3 is 
determined as g/3 = Raiisa e /(ATL 3 ). In simulations, the temperature difference AT is fixed at 10, the 
reference temperature Tq is assigned to be T c , Pr is set to be 7, Ra = 0, and Rai = 6.4 x 10 5 . In this test, 


the dimensionless parameters A and B are determined respectively by A = t/ — 0.5 — Ma^/3Pr / RaiL / 6 X 


and B = tt — 0.5 — Ma^/3/(PrRai)L/(crd x ), and the square cavity are covered by a uniform lattice of 
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Table 2: Comparisons of the average Nusselt numbers for various Ra, e and Da with Pr = 1.0 and Raj = 0 (grid sizes: 


120 X 120 for Da 

= 10 -2 

, 200 X 200 for Da = 10"' 

l , and 250 x 250 for Da = 

to- 6 ). 





Da{N x X N y ) 

Ra 


e = 0.4 



e = 0.6 



e = 0.9 


Ref. [48] 

Ref. [30] 

Present 

Ref. [48] 

Ref. [30] 

Present 

Ref. [48] 

Ref. [30] 

Present 

10" 2 

10 3 

1.010 

1.007 

1.008 

1.015 

1.012 

1.012 

1.023 

1.017 

1.018 

(120 x 120) 

10 4 

1.408 

1.362 

1.365 

1.530 

1.494 

1.498 

1.640 

1.628 

1.641 


10 5 

2.983 

3.009 

3.012 

3.555 

3.460 

3.463 

3.910 

3.939 

3.946 

10" 4 

10 5 

1.067 

1.067 

1.067 

1.071 

1.069 

1.069 

1.072 

1.073 

1.073 

(200 x 200) 

10 6 

2.550 

2.630 

2.618 

2.725 

2.733 

2.734 

2.740 

2.796 

2.818 


10 7 

7.810 

7.808 

7.811 

8.183 

8.457 

8.506 

9.202 

9.352 

9.351 

lO" 6 

10 7 

1.079 

1.085 

1.089 

1.079 

1.089 

1.094 

1.08 

1.090 

1.102 

(250 x 250) 

10 8 

2.970 

2.949 

3.014 

2.997 

2.957 

3.035 

3.00 

3.050 

3.068 


10 9 

11.460 

11.610 

11.733 

11.790 

12.092 

12.149 

12.02 

12.341 

12.399 



Figure 10: Schematic diagram of a porous cavity with internal heat-generating fluid. 


N x x N v = 120 x 120. 


In Fig. 1111 the streamlines and isotherms predicted by the present model with the Brinkman-extended 
Darcy model (e = 1 ,F e = 0) are shown for different Da. We note that these presented results are qualita- 

, . To make a quantitative comparison, the 


30 


tively consistent with those in previous numerical studies 
maximum dimensionless stream function il> m ax (normalized by Ly/gf3ATL) and the maximum dimensionless 
temperature 6 max of 9 = (T — T c )jAT are computed and listed in Table [3] From the comparisons shown in 
the Table, it is clear that our numerical results are in excellent agreement with those reported data in Refs. 
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(a) Da = oo 




(b) Da = IQ’ 2 




Figure 11: Streamlines (left) and isotherms (right) for Ra = 0 ,Raj = 6.4 X 10 s , £ = 1.0. F- = 0, and Pr = 1.0: (a) Da = oo; 
(b) Da = 10 -2 ; (c) Da = 10 -4 (gird size: 120 X 120). 
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Table 3: Comparisons of ipmax and Omax between the present results (grid size: 120 x 120) and those published solutions 
0,0 at Ra = 0 ,R ai = 6.4 x 10 5 , e = 1.0, F e = 0 , Pr = 1.0. 


Da(N x x N y = 120 x 120) 


Ref. [52] 

Ref. [30] 

Present 

oo 

'ifimax 

2.91 x 10" 3 

2.86 x 10" 3 

2.83 x 10 -3 


@max 

4.75 x 10~ 2 

4.79 x 10" 2 

4.79 x 10" 2 

10 -2 

'iftmax 

2.21 x 10" 3 

2.17 x 10" 3 

2.14 x 10~ 3 


@max 

5.22 x 10" 2 

5.26 x 10~ 2 

5.26 x 10~ 2 

10 -4 

'iftmax 

1.10 x 10 -4 

1.06 x 10~ 4 

1.04 x 10 -4 


@max 

7.35 x 10" 2 

7.34 x 10~ 2 

7.35 x 10~ 2 


3C 




52|. In addition, we can observe that as Da decreases the maximum dimensionless stream function 
c decreases, while the maximum dimensionless temperature 9 max increases. This different trend can 


be attributed to the fact that as Da decreases the flow is resisted increasingly by the presence of porous 
medium and hence the fluid velocity in the cavity decreases, while the buoyancy force is retarded to decrease 


the convective heat transfer 


521 ]. which is also indicated in Fig. [TT] 


4-4- Thermal convection in a porous cavity with internal heat generation 

In this subsection, the thermal convection flow in a porous cavity with internal heat generation and 
external sidewall heating is studied. The flow geometry and boundary conditions are identical to those 
shown in Fig. [HJ but a volumetric internal heat source Q is imposed over the domain. The temperature 
difference is AT = Th — T c , the reference temperature is To = (Tf, + T c )/2, and the internal heat source Q 
is determined by Q = Raia e AT/(RaL 2 ). In the following simulations, the Prandtl number Pr is set to be 
0.7, and the lattice size of 150 x 150 and 200 x 200 are employed for Ra = 10 5 and Ra = 10 6 , respectively. 
The parameters A and B are determined according to Eq. (15511 . and the average Nusselt number Nu are 
computed by Eq. (Rn?j) . 

Fig. ri~2l shows the streamlines and isotherms for different Rai at Ra = 10 5 ,e = 0.4, and Da = 10 -2 . At 
Rai = 10 3 , the flow fields exhibit an elliptic vortex feature in the core of the cavity, and the isotherms are 
horizontal in the cavity center due to the dominated convection effect from the external sidewall-heating. 
As Rai increases to Rai = 10 5 , the flow fields are in a similar status as those at Rai = 10 3 , but the flow 
circulation becomes more stronger. When Rai is further increased to 10 7 , due to the increase of internal¬ 
heating effect, two vortices with counter rotations exist in the cavity, and the isotherms in the left region of 
the cavity exhibit an opposite curve to those in Figs. [121(a) and (b). 

Next, the streamlines and isotherms for different Rai at Ra = 10 5 , e = 0.4 but Da = 10 -4 are depicted in 
Fig. [13] Compared with the results of Da = 10 -2 in Fig. [T2J the lower permeability of the porous medium 
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Figure 12: Streamlines (left) and isotherms (right) for Ra = 10 s , e = 0.4, Da = 10 2 , and Pr = 0.7 (gird size: 150 X 150). 
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(b) R ai = 10 5 




(c) Raj = 10' 


Figure 13: Streamlines (left) and isotherms (right) for Ra = 10 5 , e = 0.4, Da = 10 4 , and Pr = 0.7 (gird size: 


150 x 150). 
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Table 4: Comparisons of Nu and 6max between the present results and the literature results in Refs. [30|, I53l1 at Pr = 0.7 
(e = 0.4, grid sizes: 150 X 150 for Ra = 10 s , and 200 X 200 for Ra = 10 6 ). 


Ra 

Rai 


Da = 10 8 , 

e = 0.9999 

Da = 

10 -2 

Da = 

10 -4 

Da = 

10" 6 

Nu 

@max 

Nu 

@max 

Nu 

@max 

mi 

&max 

10 5 

10 3 

Ref. [53] 

4.505 

0.5 

2.884 

0.5 

1.058 

0.5 

0.995 

0.5 



Ref. [30] 

4.538 

0.5 

2.903 

0.5 

1.061 

0.5 

0.995 

0.5 



Present 

4.536 

0.5 

2.900 

0.5 

1.060 

0.5 

0.995 

0.5 


10 5 

Ref. [53] 

4.021 

0.5 

2.401 

0.5 

0.571 

0.5 

0.506 

0.5 



Ref. [30] 

4.057 

0.5 

2.421 

0.5 

0.575 

0.5 

0.506 

0.5 



Present 

4.056 

0.5 

2.421 

0.5 

0.574 

0.5 

0.506 

0.5 


10 7 

Ref. [53] 

-42.45 

5.54 

-44.08 

7.30 

-46.37 

10.86 

-48.35 

12.63 



Ref. [30] 

-41.35 

5.52 

-43.41 

7.23 

-45.08 

10.76 

-48.31 

13.27 



Present 

-40.99 

5.34 

-43.05 

7.11 

-44.79 

10.53 

-48.29 

12.62 

10 6 

10 7 

Ref. [53] 

4.129 

0.66 

1.242 

0.81 

-1.789 

1.01 

-3.93 

1.32 



Ref. [30] 

4.254 

0.65 

1.296 

0.80 

-1.728 

1.05 

-3.926 

1.32 



Present 

4.244 

0.64 

1.314 

0.80 

-1.729 

1.05 

-3.941 

1.32 


reduces the strength of the flow field here. At Rai = 10 3 , the isotherms show a weak-convection structure, 
and the center vortex with weaker circulation can be observed from the streamlines. As Raj increases to 
10 7 , the flow field is dominated by changing from the external sidewall-heating to the internal-heating, and 
a convection temperature structure results from the larger internal-heating effect. 

The effect of Da on the flow patterns and isotherms at Ra = 10 6 , Raj = 10 7 and e = 0.4 is illustrated in 
Fig. [TU It is noted for Rai/Ra > 10 that the heat transfer in the cavity is induced by the internal-heating 
more intensively than that by the external sidewall-heating. Owing to the combined influence of the internal 
heat generation and the buoyancy force, a couple of vortices appear respectively in the vicinity of the hot 
wall and the cold wall. It can be observed from the isotherms that as Da decreases, the convective heat 
transfer becomes weaker due to the lower strength of flow circulation. In addition, in the internal-heating 
dominated condition, the maximum dimensionless temperature 9 max increases with the decrease of Da. All 
the observations from Figs. H2lfl4l are in good agreement with those results in previous studies 

To quantify the results, the maximum dimensionless temperature 9 max and the average Nusselt numbers 
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a. 


along the hot wall are computed. Table 
results using the finite element method 


lists the present numerical results together with some previous 


and the LBM 130] for comparison. As can be seen from the 


Table, the present results agree quantitatively well with those reported results in the literature. 
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(b) Da = 10" 4 





Figure 14: Streamlines (left) and isotherms (right) for Ra = 10®, e = 0.4, Raj = 10 7 , and Pr = 0.7 (gird size: 200 X 200). 
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5. Conclusions 


In this paper, a modified LBGK model is proposed for simulating the incompressible fluid flow and 
convective heat transfer in porous media at the REV scale. Different from previous LBGK models for 
thermal flows in porous media, the macroscopic temperature equation can be correctly recovered from 
the proposed LBGK model by employing a modified EDF and a source term in the evolution equation of 
temperature field. In addition, following the idea of the LKS, the EDFs for the flow and temperature fields 
in the present LBGK model are incorporated with the shear rate and temperature gradient, for which a local 
computing scheme is proposed instead of the traditionally finite-difference schemes. The key advantage of 
this treatment is that the fluid viscosity and the thermal diffusivity can be determined by two additional 
parameters independent from the relaxation times, and thus better numerical stability can be achieved at 
low viscosities and thermal diffusivities. 

The modified LBGK model has been well validated by simulating four two-dimensional convective flow 
problems, including mixed heat convection in a porous channel, natural convection in a porous cavity, and 
thermal convection in porous cavity with internal heat generation with isothermally cooled walls/external 
sidewall-heating. As for the present model for the flow and temperature fields and the local schemes for 
the gradient operators, they have been demonstrated to be second-order accurate in space by the first test 
problem where the analytical solution exists. Furthermore, the computational accuracy of the present model 
and the local scheme for gradient calculations have been investigated with different values of dimensionless 
relaxation times and the two additional parameters. The results suggest that to obtain more accurate 
results, the relaxation times should be around unity, and the additional parameter in the flow field should 
be close to zero, while the other additional parameter in the temperature filed is assigned with smaller 
value as possible. In addition, the present modified LBGK model is confirmed to be more stable than the 
non-modified LBGK models at low viscosities and thermal diffusivities. Thus, the present work provides a 
more useful LBGK model in studying heat (or mass) transfer processes in porous media. 
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Appendix A. Chapman-Enskog analysis on the LBGK model for the flow field 


The Chapman-Enskog analysis is provided for the LBE (0 to recover the hydrodynamic equations. To 
this end, the following expansions with a expansion parameter £ are introduced: 


Si = if ) +£if ) +£ 2 ./f ) + ---, 

(A.la) 

9t = + £ 2 <9t 2 , 

(A.lb) 

v = £Vi, F = £F (1) . 

(A.lc) 


Note that the shear rate S resided in the EDF 0 is related with the spatial gradient. Thus, the EDF is 
expanded as the following multiscale form: 




where /, c ^ and are given by 


/■ (0) = 


Po - (1 - w 0 )fjr + PoSo(m), i = 0 

+ PO s i( M )i * 7^ 0 


/■ e (!) : (cjCj c 2 s I) 

J i &HPo 2 c 2 ’ 


(A.2) 


(A.3) 


where Si is defined as S = £Sq. With these definitions, one can easily obtain that: 

Y fi i0) = P°’ Y c i/r (0) = Po«, Y = £pI + > Y c i c i c ifi { ° ] = c sPoA • u, (A.4) 

i i i i 

Y fi {1) = °> Y c ^ {1) = °’ Y c i Ci fi {1) = clpoAdtSi, (A.5) 

i i i 

where A • u = u a 8p 1 + up5 ai + Ury5 a p, and 6 a p denotes the Kronecker tensor. 

Expanding fi(x + CiSt,t + 5t) in Eq. 0 with the Taylor theorem about space x and time i, and applying 
the above multiscale expansions to the resulting equations, one can obtain the following set of successive 
equations in the order of £: 


£° : 

AO) _ ,e(0) 

J i J i j 

(A.6a) 

£ X : 

99ii/i 0) =-^-[/f ) -/r (1) ] 

T fO t 

(A.6b) 

£ 2 : 

F) f(^) i n fW _i_ ^ x r)2 AO) _ 1 p(2) 

VtiJi +L>1 iJi • 

(A.6c) 

Substituting from Ea. (|A.6bl) into Ea. (lA.6cf) vields 


f + 1 

( ; - 2rf) + A uffl> + = -A 2K 

(A.7) 


Based on the mass and momentum conservation laws in combination with the definitions of fluid density 
and velocity, the following equations can be established: 


(A.8a) 
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(A.8b) 


p 0 u = Y c z fj; eq) = Y + ijPo F - 

% i 

Along with these results, we can get: 

£/.“’ = o <* > o), 


EE^‘ , = ° (*>n 


(A.9a) 
(A.9b) 


From Eq. for the forcing term Fi , the following velocity moments of Fi can be obtained 


Ef (1) = o, E^ 11 ’ = ( 1 - 2^ )' , ° F “ E<w*f« = (i- — 


1 \ fp 0 uF (1) poF^i 


2r fJ \ e 


(A.10) 

By taking the velocity moments of Eq. (IA.6bl) and combining Eqs. (IA.6al) . (IA.4I) . (IA.5I) . (IA.9I) and 
(IA.10II . the first-order macroscopic equations can be obtained 


5 tl p 0 + Vi- (p 0 u) = 0, 
d tl {p 0 u) + Vi • (epl + = p 0 F {1) . 

Similarly, taking the moments of Eq. (IA.7I) leads to the following equations at 0(£ 2 ): 

dt 2 Pa = 0, 

d t2 {pou) + (l- Vi • ( Y c i c *fz 1] ) + F ~' V i • (clpoAStS^ 


2 T f. 

St 

2 


1 - 


2 Tf 


Vi 


/ 2r f 

/ p 0 uFA) poF^u 

V e 


(A.11a) 
(A.lib) 

(A.12a) 

(A.12b) 


= 0. 


To proceed further, the momentum flux c i c i f- 1 ' 1 needs to be evaluated. By making use of Eqs. (IA.6bl) 
and (IA.11I) . and with some standard algebraic manipulations, we obtain that 


£ W< (1) = -TfSt Y Cid (Duf^ F ^) + Y 


= ~ T fSt 


'ti (epl + + c 2 Vi • (A • u) - fl - 


1 \ fp 0 uF (1) p 0 F^ 


u 


2 r f ) \ £ 


c 2 s p 0 AS t S i 
(A.13) 


|, we have neglected the terms of 0(Ma 2 ) from the fact that 

(A. 14) 


In the above derivations, as used in Ref. 


dp 


^ = 0(Ma 2 ), 0{Sp) = 0(Ma 2 ), 0{u) = 0{Ma). 

Here Ma is the Mach number of the flow. Therefore, the final form of Eq. (lA.12bl) can be written as 


dt 2 (Po«) + V i 


p 0 c s ( A - T f + - ) (5 t Si 


= 0. 


(A.15) 
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Combining the first- and second-order equations and noting that po is a constant, the final hydrodynamic 
equations for the flows in porous media can be derived: 


V • u = 0, 

du , „/u\ 1 „ 

— + (u • V) - =-V(ep) + ^ e v 2 ix + F, 

at \£/ po 

where v e is the effect viscosity given by 



(A.16a) 
(A.16b) 


(A.17) 


Now, we derive the formula to compute the pressure p via the distribution function. Note that vq(u ) = 0 
due to V ■ u = 0. From the expression of f^ and the Taylor expansion for the LBE (|fj]l . we then have 


t)=po- f ( Q eq) (x, t ) + p 0 S 0 (u(x, t)) 

C S 

= Po - [fo(x,t) + TfSt(Dofo(x,t) - F 0 (x,t))} + p 0 s 0 (u(x,t)) + 0(5?) 
= Po - fo(x,t) - Tf6 t d t fo{x,t) + Tf5 t F 0 (x,t) + p 0 s 0 (u(x,t)) + 0(5?). 


It can be demonstrated that the term of d t fo(x 1 1 ) in the above equation can be neglected. As for this result, 
we come back again to Eqs. (IA.2I) and (IA.3I) for f 0 (x,t), and recur to Eq. (IA.14I) to derive 


d t fo(x,t) = d t fo {0) (x,t) + 0(0 

£(l-uo)dp u} 0 p 0 d\u\ 2 
c 2 s dt 2 eel dt 

= 0(£ + Ma 2 ), 


Therefore, we can obtain 
e(l - wq) 


p(x, t)= po- f 0 (x, t ) + T f 5 t F 0 (x, t) + p 0 s 0 (u(x, t)) + 0(5? + £5 t + Ma 2 5 t ) 
= ^2 + T f 5 t F 0 (x,t) + p 0 s 0 (u(x,t)) + 0(5? +£5 t + Ma 2 5 t ), 

i^O 


(A.18) 


where Eq. (IA.8al) has been used. As a consequence, the pressure p can be calculated with accuracy of order 
0(5? + £5t + Ma 2 5t) as 


P(x,t ) = 


e(l - w 0 ) 

where sq(u) and Fq can be written in a definite form as 


5>0M) + Pos 0 (u(x,t)) + T f 5 t F 0 (x,t) 

i =£0 


(A.19) 


/ \ ^0 | |2 


. 1 \ u ■ F 

Fo = ^ o{1 - \ 


(A.20) 
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